function surface = lpcvd_2DSurface(spectrum)
fprintf(1,'Constructing lpcvd_2DSurface from spectrum2D ......\n');
h  = zeros(spectrum.LatticeSize,spectrum.LatticeSize);
dx = spectrum.XYmax/spectrum.LatticeSize;
x  = 0:dx:(spectrum.XYmax-dx);
y  = x;
mode = spectrum.mode;
for i = 1:length(x)
    fprintf(1,'Drawing surface with x = %.3f\n',x(i));
    for j = 1:length(y)
        for m = 0:mode
            for n = 0:mode
                h(i,j) = h(i,j)...
                    + phi1(spectrum,m,n,x(i),y(j))...
                    + phi2(spectrum,m,n,x(i),y(j))...
                    + phi3(spectrum,m,n,x(i),y(j))...
                    + phi4(spectrum,m,n,x(i),y(j));
            end
        end
    end
end
surface = lpcvd_2DSurface(h,x,y,'nm');